#include "plant2D_EW.h"

const double k_w = 3.3829e-12;
const double a_w = 0.6042;
const double E_w = 2.7e-3; // Unit = EV
const double k_c = 1.0274e-13;
const double a_c = 0.1669;
const double E_c = 1.9e-3; // Unit = EV
const double a_v = 15.55493;
const double k_v = 20.64504;
const double a_t = 0.02332;
const double k_t = 0.0261;
const double k_B = 8.617343e-5; // Boltzmann's constant, Unit = ev/K

// --------------------------------------------------------------
plant_SPDE2D::plant_SPDE2D(int M, int L,double idt){
  int i,m,n;
  mode = M;	
  int idim[2] = {1,1}; // Input#1 = T(K), Input#2 = W (layer/s)
  int odim[2];
  odim[0] = (mode+1); // The mode with sin*sin
  odim[1] = (mode+1); // 

  block::initialize(2,idim,2,odim,idt);

  LatticeSize = L;

  z1 = block::outputs_next[0];
  z2 = block::outputs_next[1];

  lambda = new double*[mode+1];
  for(m=0;m<mode+1;m++){
    lambda[m] = new double[mode+1];
  }

  for(i=0;i<2;i++){
    Km2[i] = new double[(mode+1)];
  }

  double dphi;
  int ix,ix_next;
  int index;
  for(m=0;m<=mode;m++){
 
      index = index1D(m);
      for(i=0;i<2;i++){
	Km2[i][index] = 0.0;
      }
      for(ix=0;ix<LatticeSize;ix++){
	  if(ix == LatticeSize-1){ix_next = 0;}
	  else{ix_next = ix+1;}
	  dphi = phi1m(m,ix_next)-phi1m(m,ix);
	  Km2[0][index] += dphi*dphi;
	  dphi = phi2m(m,ix_next)-phi2m(m,ix);
	  Km2[1][index] += dphi*dphi;
      }
      Km2[0][index] = Km2[0][index]/(pi*pi);
      Km2[1][index] = Km2[1][index]/(pi*pi);
    }
  }
  /*	ofstream fout_Km2;
  fout_Km2.open("Output_Km2.dat");
  for(m=0;m<=mode;m++){
  for(n=0;n<=mode;n++){
  fout_Km2 << Km2[1][index2D(m,n)] << " ";
  }
  fout_Km2 << "\n";
  }
  fout_Km2.close();*/

  dt_con = 0.0001;
  sqrt_dt_con = sqrt(dt_con);
  dt_conservative_factor = 0.1;

  gsl_rng_env_setup();
  noise_type = gsl_rng_default;
  noise = gsl_rng_alloc(noise_type);

  reset();
}
// --------------------------------------------------------------
plant_SPDE2D::~plant_SPDE2D(){
  int m;
  for(int i=0;i<4;i++)
  {
    delete []Km2[i];
  }

  for(m=0;m<mode+1;m++){
    delete []lambda[m];
  }
  delete []lambda;
  gsl_rng_free(noise);
}
// --------------------------------------------------------------
void plant_SPDE2D::reset(){
  block::time = 0.0;
  int m,n;
  for(m=0;m<mode+1;m++){
    for(n=0;n<mode+1;n++){
      z1[index2D(m,n)] = 0.0;
      z2[index2D(m,n)] = 0.0;
      z3[index2D(m,n)] = 0.0;
      z4[index2D(m,n)] = 0.0;
    }
  }

  T = 610;
  W = 0.5;
  c = eq_c(T,W);
  c2 =eq_c2(T,W);
  sigma = eq_sigma2(T,W);

  for(m=0;m<mode+1;m++){
    for(n=0;n<mode+1;n++){
      lambda[m][n] = -4*c2*(m*m+n*n);
    }
  }
}
// --------------------------------------------------------------
void plant_SPDE2D::update(double sysTime){
  double T = inputs[0][0];
  double W = inputs[1][0];
  update(sysTime,T,W);
}
// --------------------------------------------------------------
void plant_SPDE2D::update(double sysTime, double iT, double iW){	
  if(time>sysTime) return;

  int m,n;
  int index;
  double c_200 = c/(LatticeSize*LatticeSize);

  while(time<=sysTime){
    if(iT!=T || iW!=W){
      T = iT;
      W = iW;
      c = eq_c(T,W);
      c2 = eq_c2(T,W);
      sigma = sqrt(eq_sigma2(T,W));
      for(m=0;m<mode+1;m++){
	for(n=0;n<mode+1;n++){
	  lambda[m][n] = -4*c2*(m*m+n*n);
	  dt_con[m][n] = -2/lambda[m][n]*dt_conservative_factor;
	}
      }

    }
    for(m=0;m<=mode;m++){
      for(n=0;n<=mode;n++){
	index = index2D(m,n);
	if(m>0&&n>0){
	  z1[index] = z1[index]+lambda[m][n]*z1[index]*dt_con+gsl_ran_gaussian_ziggurat(noise,1.0)*sigma*sqrt_dt_con;
	}

	if(m==0&&n==0){
	  z2[index] = z2[index]+c_200*dt_con+gsl_ran_gaussian_ziggurat(noise,1.0)*sigma*sqrt_dt_con;
	}
	else{
	  z2[index] = z2[index]+lambda[m][n]*z2[index]*dt_con+gsl_ran_gaussian_ziggurat(noise,1.0)*sigma*sqrt_dt_con;
	}

	if(m>0){
	  z3[index] = z3[index]+lambda[m][n]*z3[index]*dt_con+gsl_ran_gaussian_ziggurat(noise,1.0)*sigma*sqrt_dt_con;
	}

	if(n>0){
	  z4[index] = z4[index]+lambda[m][n]*z4[index2D(m,n)]*dt_con+gsl_ran_gaussian_ziggurat(noise,1.0)*sigma*sqrt_dt_con;
	}
      }

    }
    time += dt_con;
  }
  block::output();
  return;
}
// --------------------------------------------------------------
void plant_SPDE2D::output(string FileName){
  ofstream fout;	
  fout.open(FileName.c_str());
  fout << "Mode = " << mode;
  int i;

  int dim = (mode+1)*(mode+1);

  for(i=0;i<dim;i++){
    fout << z1[i] << " ";
  }
  fout << "\n";

  for(i=0;i<dim;i++){
    fout << z2[i] << " ";
  }
  fout << "\n";

  for(i=0;i<dim;i++){
    fout << z3[i] << " ";
  }
  fout << "\n";

  for(i=0;i<dim;i++){
    fout << z4[i] << " ";
  }
}
// --------------------------------------------------------------
double plant_SPDE2D::getM2(){
  int m,n,index;
  double m2 = 0.0;
  for(m=0;m<=mode;m++){
    for(n=0;n<=mode;n++){
      index = index2D(m,n);
      m2 +=  Km2[0][index]*z1[index]*z1[index]
      +Km2[1][index]*z2[index]*z2[index]
      +Km2[2][index]*z3[index]*z3[index]
      +Km2[3][index]*z4[index]*z4[index];
    }
  }
  return m2;
}
// --------------------------------------------------------------
double plant_SPDE2D::getR2(){
  int m,n,index;
  double r2 = 0.0;
  for(m=0;m<=mode;m++){
    for(n=0;n<=mode;n++){
      if(m!=0||n!=0){
	index = index2D(m,n);
	r2 +=  z1[index]*z1[index]
	+z2[index]*z2[index]
	+z3[index]*z3[index]
	+z4[index]*z4[index];
      }
    }
  }
  r2 /= (LatticeSize*LatticeSize);
  return r2;
}
// --------------------------------------------------------------
double plant_SPDE2D::eq_c(double T,double W)
{
  return W*(1-k_w*pow(W,-a_w)*exp(k_B*T/E_w));
}

double plant_SPDE2D::eq_c2(double T,double W)
{
  return k_c*exp(k_B*T/E_c)*pow(W,-a_c)/(LatticeSize*LatticeSize);
}

double plant_SPDE2D::eq_sigma2(double T,double W)
{
  double index = (a_t+k_t*W)*T-a_v-k_v*W;
  return pi*pi/(LatticeSize*LatticeSize)*(1+exp(index));
}
